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Section  I 


INTRODUCTION 

The  USA-STAG S  code  [1]  computes  the  transient  response  of 
a  submerged  structure  that  is  subject  to  underwater  shock  ex¬ 
citation.  The  structural  behavior  may  be  linear  or  nonlinear. 
The  success  of  this  code  in  treating  the  complex  phenomena  in¬ 
volved  is  largely  due  to  the  development  of  the  Doubly  Asymptotic 
Approximations  [2,3]  that  describe  the  fluid-structure  interaction 
in  terms  of  variables  at  the  “wet”  surface  of  the  structure  only, 
thus  avoiding  the  modeling  and  computational  burden  produced 
by  surrounding  the  structure  with  fluid-volume  elements.  This 
approach  has  opened  the  way  to  greater  understanding  of  the 
problem  area  through  computer  simulation.  A  restricted  class  of 
problems  exist,  however,  for  which  the  possible  effects  of  hull 
cavitation  on  the  structural  response  must  be  considered.  A 
related  question  concerns  the  treatment  of  internal  fluids  con¬ 
tained  within  ballast  tanks,  free-flooded  areas,  etc.  An  exten¬ 
sion  of  the  present  underwater-shock  computational  technology 
to  include  fluid- volume  elements  is  required  to  treat  both  of  these 
complications. 

A  highly  efficient  computational  scheme  for  treating  a  cavitat- 
ing  acoustic  fluid  has  been  devised  by  Newton  [4].  The  scheme 
involves  the  use  of  the  displacement  potential,  which  is  a  scalar 
quantity,  as  the  primary  variable  in  the  formulation  of  finite- 
element  matrix  equations  for  the  fluid  volume.  This  choice  has 
significant  advantages  over  a  displacement  vector  formulation, 
which  triples  the  number  of  fluid- volume  unknowns  and  does  not 
automatically  enforce  irrotationality  of  fluid  motions.  A  second 
feature  of  Newton’s  scheme  is  the  use  of  an  explicit  time  integra- 


tion  method.  These  two  concepts  form  the  basis  of  the  newly 
constructed  Cavitating  Fluid  Analyzer  (CFA)  that  has  been  in¬ 
terfaced  with  USA-STAGS. 

The  theory  behind  the  implementation  of  CFA  and  its  in¬ 
teraction  with  both  USA  and  STAGS  is  presented  in  Section  n, 
along  with  stability  analyses  of  the  staggered  time-integration 
procedure  used  in  the  coupled  USA-STAGS-CFA  system.  Section 
in  contains  a  brief  discussion  of  the  software  implementation  and 
usage  of  this  system.  Finally,  Section  IV  presents  computational 
results  for  two  simple  cavitation  problems  for  which  solutions  are 
available,  namely,  the  one-dimensional  Bleich-Sandler  flat-plate 
problem  [5],  and  a  variant  of  a  two-dimensional  cylindrical-shell 
problem  studied  by  Newton  [6].  The  USA-STAGS-CFA  results  are 
shown  to  be  in  excellent  agreement  with  the  previous  solutions. 

It  should  be  emphasized  that  this  report  is  not  a  users  manual 
for  the  USA-STAGS-CFA  system.  Such  a  manual  will  be  issued 
once  the  new  software  has  undergone  further  evaluation  in  actual 
three-dimensional  problems  and  has  been  provided  with  a  CFA 
user  interface  that  meets  the  production-level  needs  of  the  under¬ 
water  shock  community. 


Section  II 


THEORY 


§2.1  Problem  Description 

A  structure  is  submerged  in  a  fluid  idealized  as  an  infinite  acous¬ 
tic  medium  incapable  of  transmitting  tensile  stresses.  A  compres¬ 
sive  shock  wave  propagates  through  the  fluid  and  impinges  on 
the  structure.  If  the  structure  is  sufficiently  flexible  and  the  am¬ 
bient  hydrostatic  pressure  sufficiently  low,  the  scattered  negative 
pressure  wave  may  induce  cavitation  in  the  subregion  that  was 
traversed  by  the  incident  shock  wave  before  reaching  the  struc¬ 
ture.  This  phenomenon  is  known  as  hull  cavitation. 

Because  of  the  nonlinear  nature  of  cavitation,  a  boundary- 
element  treatment  of  the  entire  fluid  domain  as  a  “DAA  mem¬ 
brane”  surrounding  the  structure  is  ruled  out.  (Boundary  element 
methods  are  restricted  to  homogeneous  linear  domains.)  Instead, 
a  realistic  computer  analysis  of  this  problem  requires  the  con¬ 
sideration  of  the  three  interacting  fields  illustrated  in  Figure  1: 
structure  S,  cavitating  fluid  volume  V,  and  DAA  membrane  D. 

The  DAA  boundary  should  be  placed  as  far  away  as  neces¬ 
sary  to  encompass  the  cavitating  fluid  subregion.  Inasmuch  as 
the  extent  of  the  latter  is  generally  unknown  before  the  analysis 
is  performed,  some  iterations  on  the  placement  of  the  DAA  boun¬ 
dary  may  be  necessary.  For  example,  if  an  analysis  shows  cavita¬ 
tion  occurring  at  points  on  the  DAA  boundary,  the  latter  should 
be  moved  further  away  from  the  structure  interface.  Conversely, 
if  a  preliminary  analysis  shows  that  wide  external  fluid  regions 
remain  pressurized,  the  DAA  boundary  might  be  moved  closer  to 
the  structure  to  reduce  the  computational  cost. 
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The  nature  of  the  interaction  among  the  three  fields  (5,  V, 
D )  can  be  pictorially  illustrated  as  follows. 


displacements 

pressures 

=> 

=> 

Submerged 

Fluid 

DAA 

Structure 

<= 

Volume 

<= 

Boundary 

pressures 

displacements 

This  diagram  shows  that  the  fluid  volume  basically  functions 
as  a  pressure  transducer  between  the  DAA  boundary  and  the 
submerged  structure.  By  way  of  contrast,  the  more  conventional 
two-field  structure-DAA  interaction  can  be  diagrammed  as 


velocities 


Submerged  DAA 

Structure  Boundary 

<= 

pressures 

This  section  presents  theoretical  background  material  with 
emphasis  on  the  finite  element  discretization  of  the  fluid-volume 
subsystem  and  the  treatment  of  the  coupled  problem  with  a  stag¬ 
gered  time-integration  procedure. 
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§2.2  Fluid- Volume  Field  Equations 


This  subsection  is  a  compendium  of  well  known  continuous  field 
equations  for  linear  and  bilinear  acoustic  fluids.  These  equations 
are  collected  here  to  introduce  notation  and  to  make  this  report 
reasonably  self-contained. 

Small,  Irrotatioual  Motions  of  a  Compressible  Fluid.  Let  pH 
and  p  be  the  hydrostatic  pressure  and  fluid  density,  respectively, 
in  the  acoustic  fluid  that  occupies  the  volume  V  of  Figure  1. 
Compressive  pressures  are  conventionally  denoted  as  positive.  Ma¬ 
terial  points  in  V  are  identified  by  the  global-coordinate  vector 
£  —  ( X.Y ,Z ).  If  /  =  J(X)  is  the  body  force  field,  then  the 
static  equilibrium  vector  equation  is 


Vpf/  +  ?  =  3.  (2.1) 

Let  a!  =  2(^)  be  the  fluid-particle  displacement  field  under 
dynamic  conditions.  Then 

3  =  2-2h,  (2.2) 

is  the  fluid  particle  displacement  relative  to  a  reference  hydro¬ 
static  displacement  .  The  dynamic  vector  equation  of  motion 
is 

p2  =  ip-\-J,  (2.3) 

where  p  is  the  total  pressure. 

For  small  irrotational  fluid  motions,  the  field  2  is  derivable 
from  the  scalar  displacement  potential  $  defined  by 

V#  =  ~p2.  (2.4) 

The  factor  p  is  introduced  in  the  potential  definition  for  notational 
convenience. 
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Tiuffl  tz.l/,  \z.o)  SOT  (2.4)1 

V#  =  V(p  —  pH).  (2.5) 

which  can  be  spatially  integrated  to  yield 

*  =  V-VH,  (2-6) 

For  future  use,  define  the  “densified”  relative  condensation  as 

5  =  -pvl  (2.7) 

From  (2.4)  and  (2.7)  it  follows  that 

s  =  V2^.  (2.8) 


Linear  Fluid.  The  constitutive  equation  of  a  linear  acoustic  fluid, 
valid  for  s  <  p,  is 

p—pn—c2s.  (2.9) 

Here,  c  is  the  reference  sound  speed,  which  is  linked  to  the  bulk 
modulus  K  and  density  p  by  the  relation  c2  —  K Ip. 

Bilinear  Fluid.  In  real  fluids,  cavitation  is  a  microscopically  hetero¬ 
genous  phenomenon  influenced  by  gas  dilution  concentrations.  A 
simple  yet  effective  mathematical  model  for  this  phenomenon  con¬ 
sists  of  assuming  that  the  cavitating  region  is  macro scopically 
homogeneous  and  at  zero  total  pressure.  This  leads  to  the  notion 
of  a  bilinear  fluid,  whose  constitutive  properties  are  adjusted  so 
that  it  cannot  transmit  negative  total  pressures.  The  constitutive 
equation  of  the  bilinear  fluid  is 


if  5  >  — pH l c 2, 


p  =  0 


otherwise. 


(2.10) 


The  equation  of  motion  (2.6)  becomes 


&  =  p  —  pH  if  s>—pH/c2, 

&  =  —pH  otherwise.  (2.11) 

Remark  /.  Temporal  differentiation  of  (2.4)  furnishes  the  velocity  potential 
equation 

=  —pH,  (2.12) 

« 

where  3  =  3  is  the  fluid-particle  velocity.  However,  the  velocity  potential  is 
not  a  primary  field  variable  in  the  present  work.  The  use  of  as  primary 
variable  has  been  found  (4,7)  to  produce  temporal  discontinuities  that  are 
avoided  by  the  displacement  potential  formulation. 

Remark  2.  Elimination  of  p  and  s  from  (2.6),  (2.8)  and  (2.9)  shows  that,  for 
a  linear  fluid,  ^  satisfies  the  scalar  wave  equation 

tfr  =  c2V2V  (2.13) 

and,  of  course,  so  do  p  and  s. 
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§2.3  Fluid- Volume  Discretization 


Acoustic  Finite  Elements .  The  volume  V  occupied  by  the  cavitat- 
ing  fluid  is  divided  into  three-dimensional  finite  element  sub¬ 
domains  interconnected  at  element  node  points.  The  basic  ele¬ 
ment  used  in  this  study  is  the  eight- point  isoparametric  “brick”, 
which  is  a  hexahedron  whose  geometry  is  uniquely  defined  by  the 
position  of  its  eight  corner  points.  Six-node  “wedge”  elements 
may  also  be  used  for  things  like  rounding  corners,  etc. 

Following  standard  finite  element  techniques,  the  geometry  of 
the  discretized  fluid  volume  is  expressed  in  the  standard  matrix 
form 

X  =  N‘X,  Y  =  N*  Y,  Z  =  N‘Z,  (2.14) 

where  X,  Y  and  Z  are  column  vectors  of  nodal  values  of  the 

global  coordinates  =  (X,Y,Z),  N  is  a  column  vector  of 
finite  element  shape  functions  associated  with  these  nodes,  and 
superscript  t  denotes  transposition. 

Matrix  Equations.  The  displacement-potential  Galerkin  formula¬ 
tion  of  R.  E.  Newton  [4,7]  is  used  to  derive  the  finite  element 
matrix  equations.  The  primary  field  variables  of  this  formulation 
are  the  displacement  potential  &  and  the  condensation  s.  In  ac¬ 
cordance  with  the  isoparametric  finite  element  concept,  these  two 
fields  are  interpolated  with  the  geometry  shape  functions: 


<P  =  N‘  s  =  N‘  8,  (2.15) 

where  and  s  are  the  column  vectors  of  node  values  of  &  and  5, 
respectively. 
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The  discrete  counterpart  of  the  strain-displacement  equation 
(2.8),  i.e. 

s  -  V2!?  =  0,  (2.16) 


is  the  Galerkin  equation 


/, 


N(s-  V2<P)dV  =  0. 


(2.17) 


Application  of  the  divergence  theorem  to  (2.17)  produces 

I  (Ns  +  VNV!P)dV  =  [  N ~dB,  (2.18) 

Jv  Jb  on 

where  n  denotes  the  outward  normal  to  the  boundary  B  of  V. 
Insertion  of  the  finite  element  interpolation  assumptions  (2.15) 
into  the  left-hand  side  of  (2.18)  yields  the  Galerkin  matrix  equa¬ 
tion 

Qs  =  -H*  +  b.  (2.19) 

In  this  equation  Q  and  H  are  symmetric  square  matrices  given  by 


Q=  f  NN ldV, 
Jv 


H 


=  /  (VN)(VN)ed^, 
Jv 


(2.20) 


(2.21) 


while  the  column  vector  b  is  defined  by 


=  / N 


d# 


b  dn 


dB. 


(2.22) 


The  determination  of  the  entries  of  Q  and  H  is  performed  by  stan¬ 
dard  numerical  integration  techniques  for  isoparametric  elements, 
using  Gauss-Legendre  quadrature  rules.  For  the  time-marching 
calculations  described  in  following  sections,  the  entries  of  matrix 
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H  need  never  be  explicitly  calculated;  instead,  the  vector 


r  =  — Htf  (2.23) 

is  formed  directly  in  the  numerical  integration  loop. 

The  boundary  interaction  vector  b  can  be  split  into 

b  =  b'  +  b<',  (2.24) 

where  terms  bs  and  bd  come  from  the  contributions  of  the  struc- 
ture  and  DAA  boundary,  respectively.  The  calculation  of  these 
terms  is  discussed  in  §2.5  and  §2.6. 

Dynamic  Equations.  The  space-discrete  counterpart  of  the  field 
equation  of  motion  (2.6)  is 

*  =  P~PW,  (2.25) 

where  p  and  pH  are  column  vectors  of  node  values  of  p  and  pH , 

..  - 1 -- 
»  c.  prvt  i  vnj  . 

Eqs.  (2.19)  and  (2.25)  may  be  combined  by  eliminating  p  and 
s  to  yield 

*  +  c2Q-1  H  =  c2Q_1  b.  (2.26) 

This  is  the  finite  element  discretization  of  the  wave  equation  (2.13) 
with  a  forcing  boundary  term.  (This  term  results  from  the  use  of 
the  divergence  theorem.) 

In  the  time  integration  scheme  described  in  §2.4,  the  com¬ 
bined  equation  (2.26)  is  not  used.  Instead,  it  is  far  more  con¬ 
venient  to  use  the  two  equations  (2.19)  and  (2.25)  in  tandem,  as 
the  intermediate  vector  quantities  s  and  p  are  of  interest  for  test¬ 
ing  the  cavitation  condition  (1.14)  and  for  the  determination  of 
the  interface  forces  that  act  on  the  structure  and  DAA  boundary. 
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Capacitance  Lumping.  Matrix  Q  is  the  analog  of  the  consistent 
capacitance  (mass)  matrix  of  finite  element  thermal  (mechanical) 
models.  For  explicit  time-marching  calculations  it  is  convenient 
to  replace  Q  by  a  “lumped  capacitance”  diagonal  matrix  Q  to 
avoid  solving  systems  of  linear  equations.  Thus  (2.19)  becomes 

Q  s  =  — H  #  -f-  b.  (2.27) 

Lumping  is  effected  by  assigning  the  row  sums  of  Q  to  the  diagonal 
of  Q. 


Remark  1.  Some  differences  with  Newton’s  formulation  as  presented  in  Refs. 
[4,7]  should  be  noted.  Newton  uses  a  displacement  potential  9  defined  in 
terms  of  total  displacements,  so  that  the  s  corresponding  to  (2.7)  is  the 
absolute  condensation.  Furthermore,  he  selects  the  opposite  sign  in  the 
displacement  potential  definition;  thus  his  Eqs.  (2.6)  and  (2.8)  are  &  =  pH — p 
and  s  —  — V2!?,  respectively. 

Remark  2.  Comparing  (2.26)  and  (2.13)  it  is  plain  that  — Q_1  H  is  a  dis¬ 
cretization  of  the  Laplacian  operator  V2.  and  so  is  the  “lumped  capacitance” 

form  — Q  1 H.  A  comparison  of  these  discrete  finite  element  operators  with 
conventional  finite  difference  “molecules”  is  illuminating.  At  interior  points 
of  a  regular  two-dimensional  mesh  of  square  cells,  — Q_1H  is  equivalent  to 
the  conventional  5-point  “Laplacian  star”  finite  difference  molecule.  On  the 

other  hand,  the  lumped  form  — Q  1  H  turns  out  to  correspond  to  a  9-point 
molecule  obtained  by  averaging  two  standard  5- point  molecules,  one  being 
rotated  by  45°.  A  similar  but  more  complicated  analogy  holds  for  three 
dimensions. 
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§2.4  Time  Iutegration 

For  numerical  stability  reasons  discussed  later  in  this  and  follow¬ 
ing  subsections,  the  introduction  of  artificial  damping  terms  in  the 
fluid-volume  equations  of  motion  turns  out  to  be  highly  beneficial. 
To  simplify  matters,  however,  we  start  the  presentation  of  fluid- 
volume  time  integration  techniques  assuming  that  such  terms  are 
absent. 

Undamped  Integration.  Equation  (2.25)  is  numerically  integrated 
in  time  with  an  explicit  central  difference  scheme.  Let  (.)n 
denote  the  computed  value  of  (.)  at  the  n-th  time  station  tn. 
Computations  have  proceeded  until  tn.  If  h  is  the  time  increment 
f«+ 1  —  tn,  the  central  difference  scheme  can  be  written 

'j'n+I/2  =  'l'n-l/2  +  /l(p„-pW))  (2.28) 

*»+i  =  *«  +  **„+,„.  (2.29) 

The  advancing  step  is  completed  by  the  discrete  strain-displace- 
ment  equation  (2.19)  and  the  constitutive  equation  (2.9): 

2.27 

Q*n+1  =  — 4-  bn+i,  (2.30) 

Pn+i  =  pw  +  c2in+1.  (2.31) 

Inasmuch  as  Q  is  a  diagonal  matrix,  the  solution  of  the  linear 
system  (2.30)  for  s„_j_ !  is  trivial. 

Note  that  in  this  process  s  and  p  are  computed  at  full 
stations,  whereas  4*  is  computed  at  half  stations: 

*»-  i  *n+l 

o  o 

^n  — 1/2  ^n  +  1/2  "* 

Pn  -  1  Pn  Pn-f-1 
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This  has  important  implications  in  the  starting  procedure 
described  in  §2.7. 

If  the  forcing  vector  b  does  not  depend  on  left-hand  side 
variables,  the  numerical  stability  of  the  central  difference  scheme 
(2.26)-(2.31)  is  well  known.  It  can  be  expressed  as  the  Courant- 
Friedrichs-Levy  (CFF)  condition 

h  <  hc  =  L/c ,  (2.32) 

where  L  is  the  smallest  cross  dimension  of  a  fluid-volume  finite 
element;  he  is  called  the  Courant  timestep.  The  CFL  condition 
(2.32)  will  be  derived  later  in  this  section  as  a  particular  case  of 
the  damped  system  integration. 

The  accuracy  properties  of  the  central  difference  operator  are 
also  well  known.  It  does  not  introduce  numerical  damping.  On 
a  regular  grid  of  equal-side  elements,  integrating  with  h  =  he 
furnishes  exact  nodal  results  for  the  propagation  of  plane-waves  of 
rectangular  profile  (e.g.  step-waves)  along  gridlines.  For  all  other 
cases  (h  ^  he,  irregular  grids,  or  general  waveforms),  numerical 
disoersion  occurs  and  the  solution  is  inexact. 

Artificial  Damping .  Newton’s  studies  [4,7]  have  indicated  that 
occurrence  of  cavitation  can  induce  growing  spurious  pressure 
oscillations.  These  oscillations  eventually  cause  “fragmentation” 
of  the  cavitation  region  in  the  sense  that  small  pressurized  islands 
appear  in  the  cavitation  region  while  small  zero-pressure  bubbles 
appear  in  the  pressurized  region.  The  phenomenon  has  been 
termed  frothing. 

The  most  effective  cure  to  frothing  is  numerical  damping  that 
increases  with  frequency.  This  can  be  achieved  by  augmenting 
(2.25)  with  an  artificial  damping  term  proportional  to  8: 

4-  =  p  —  pH  +  Phc2  S,  (2.33) 

in  which  (3  is  a  dimensionless  damping  coefficient  that  varies  from 
0  to  1.  The  value  of  §  is  estimated  by  a  backward  difference 
formula. 
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For  a  pressurized  fluid  (2.33)  would  reduce  to  the  slightly 
simpler  form 

^  =  P  -  P H  +  Php,  (2.34) 

since  §  =  p/c2.  But  the  expression  (2.33)  is  more  general,  as  it 
can  be  used  for  a  cavitation  region  in  which  £  is  unrelated  to  §. 
The  modified  solution-advancing  process,  including  the  test 


for  cavitation,  is  as  follows. 

—  (Sn  Sn_x)//l,  (2.35) 

^n  +  1/2  =  — 1/2  +  h  (p n  —  f>H  +  fihc2  §n),  (2.36) 

^n+1  =  “h  h  (2.37) 

Qs„-t-i  =  — H^n+i  +  bn+i,  (2.38) 

Pn+i  =  max(pH  H-  c2sn+i,0).  (2.39) 


Stability  of  Damped  Integration.  The  following  study  investigates  how  the 
stable  timestep  depends  on  the  artificial  damping  coefficient,  which  is  ob¬ 
viously  a  question  of  practical  importance.  The  effect  of  structure-  and  De¬ 
coupling  terms  on  numerical  stability  is  ignored  here.  These  two  effects  are 
examined  in  §2.5  and  §2.6,  respectively. 

The  stability  analysis  of  the  advancing  process  (2.35)-(2.39)  is  under¬ 
taken  with  a  Fourier  (normal  mode)  method,  which  involves  fairly  conven¬ 
tional  steps.  First,  the  nonhomogeneous  vector  terms:  hydrostatic  pressure 
p11  and  boundary  force  b  are  discarded  (the  latter  because  of  the  interaction 
neglect  as  noted  previously).  Next,  the  state  variables  (tf,  s  and  p  )  are 
expanded  into  normal-mode  motions  associated  with  the  eigenproblem 

Q  v  =  X  Hv  (2.40) 

Inasmuch  as  both  Q  and  H  are  symmetric  and  nonnegative  definite  (Q  is  in 
fact  positive  definite),  all  eigenvalues  X  of  (2.40)  are  real  and  positive. 
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The  advancing  step  (2.35)- (2.39)  is  rewritten  for  a  specific  normal  mode 
of  eigenfrequency  X : 

IPn  +  i/2  ~  ^n  — 1/2  +  h  IP"  +  0C2isn  ~  *n-l)j.  (2-41) 

V'n-f-l  =  V'n  +  ^$„+i/2'  (2.42) 

Sn  +  1  =  — X^n  +  1,  (2-43) 

Pn  +  1  =  lC2Sn  +  l  =  — 7*C2V'n+ 1-  (2.44) 

In  these  equations,  \p,  s  and  p  denote  amplitudes  of  ¥,  •  and  p,  respectively, 
and  7  is  a  switch  variable  that  takes  the  value  1  for  pressurized  fluid  (p  >  0) 
and  0  for  cavitating  fluid  (p  ==  0).  (Implicit  in  the  use  of  7  is  the  assumption 
that  a  normal  mode  pertains  wholly  to  either  condition;  this  can  be  only 
justified  a  posteriori  by  showing  that  only  one  condition  is  critical.) 

Elimination  of  the  three  intermediate  variables  tp,  s  and  p  yields  the 
difference  equation 

r/jn  +  l  -  Wn  +  V'r.-l  =  -d(T  +  W»  ~  Mn- 1],  (2  45) 

with 

(  =  AVX,  (2.46) 


wilicii  is  it  uimcusiuuiess  parauietei. 

The  associated  characteristic  polynomial  in  the  complex  amplification 
variable  z  (the  discrete  Laplace  transform  image  of  d>)  is 

C(z )  =  z2-l 2  -  ?(7  +  P)]z  +  1  -  (2  47) 

The  corresponding  Routh  polynomial,  obtained  through  the  involutory 
mapping  z  =  (y  - f-  1  )/{y  —  1),  is 

f (7  +  P)  y2  +  2f/?  y  +  4  -  f(7  +  2/?),  (2.48) 

from  which  it  is  easy  to  deduce  the  stability  conditions 

h  <  —  2  ■  —  and  p  >  0.  (2.49) 

c>/X(7  +  2/?) 

This  expression  shows  that  the  smallest  stable  timestep  is  associated 
with  the  largest  eigenvalue  X  of  (2.40).  It  can  be  shown  that  for  trilinear  shape 
functions  X  is  bounded  above  by  4/L2,  where  L  is  the  smallest  finite  element 
mesh  dimension.  (This  eigenvalue  is  associated  with  “hourglass"  geometric 
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distortions.)  Inserting  this  bound  in  (1.49)  yields  the  desired  formula 

h  <  hmax  =  — — - ,  (2.50) 

\A  +  20 

where  he  is  the  Courant  timestep  defined  by  (2.32). 

Setting  7  =  0  increases  the  stability  limit  if  /?  >  0.  Thus  the  occurrence 
of  cavitation  does  not  have  a  detrimental  effect  on  stability.  From  now  on 
we  can  conservatively  set  7  =  1. 

If  P  =  0,  the  stability  limit  (2.32)  of  the  undamped  central  difference 
scheme  (2.28)-(2.31)  results.  But  if  0  >  0,  the  stability  limit  is  reduced  by 

the  factor  1  / \J\  -f-  2/?.  This  factor  reaches  l/>/3  =  0.577  for  (5  =  1,  which 
is  the  maximum  suggested  damping. 

There  is  in  fact  a  slightly  smaller  recommended  maximum  timestep. 
This  is  the  transition  stepsize  htr  at  which  the  roots  of  the  quadratic  charac¬ 
teristic  equation  C{z)  =  0  pass  from  imaginary  to  real.  Imposing  the  double 
root  condition  readily  yields 


htr  — 


he 

l+P~ 


(2.51) 


This  is  smaller  than  hmax,  but  differs  by  at  most  15%  from  it  if  ()  <  1,  as 
illustrated  by  the  following  table. 


p 

hmax  f  he 

htr /he 

0.000 

1.000 

1.000 

0.250 

0.816 

0.800 

0.500 

0.707 

0.667 

0.750 

0.632 

0.571 

1.000 

0.577 

0.500 

For  plots  of  hjhc  vs.  f),  see  Figures  2  through  5  in  §2.5. 
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§2.5  Structure-Fluid  Interaction 

Action  of  Fluid-Volume  on  Structure.  The  differential  equation 
of  motion  for  the  dynamic  response  of  a  structure  spatially  dis¬ 
cretized  by  the  finite  element  method  can  be  expressed  in  the 
form 

Msr  +  Csr  +  Ksxs  =  fst  (2.52) 

where  xs  is  the  column  vector  of  nodal  structure  displacements, 
M<.,  Cs  and  Ks  are  the  structural  mass,  damping  and  stiffness 
matrices,  respectively,  and  fs  is  the  external  nodal  force  vector. 
For  acoustic-wave  excitation  of  a  submerged  structure  through 
the  fluid-volume  mesh,  f  is  given  by 


f  s  —  GgAsP,  (2.53) 

where  p  is  t  he  column  vector  of  total  pressures  at  the  nodes  of  the 
fluid- volume  mesh,  as  in  (2.25),  As  is  a  diagonal  matrix  of  con- 
IrihutiriP-  surface  nrens  fuirrnunrlina  fluid- volume  nodes  in  contact 
with  the  structure,  and  Gs  is  the  transformation  matrix  that  re¬ 
lates  structure  and  fluid  nodal  surface  forces.  If  the  structure  and 
fluid-volume  nodes  are  in  one-to-one  correspondence,  Gs  reduces 
to  the  identity  matrix  for  all  wet-surface  structure  nodes,  and  is 
zero  ot  herwise. 

Action  of  Structure  on  Fluid-Volume.  The  effect  of  the  structural 
response  on  the  fluid  volume  field  resides  in  the  boundary  inter¬ 
action  term  b*  implicitly  defined  by  (2.23)  and  (2.24): 

b9  =  /  N  —  dB,  (2.54) 

Jbs  dn 

where  is  the  wet  structure  (contact)  surface.  To  evaluate  bp, 
replace 

d* 

_  ^  -pw  =  pr  n, 


2-16 


(2.55) 


and 


xs  =  Ng  xs.  (2.56) 

in  (2.54).  Here  vj  is  the  structural  displacement  normal  to  the 
structure’s  wet  surface;  fl  is  the  wet  surface  normal  vector  con¬ 
sidered  as  positive  going  into  the  fluid,  and  Ns  is  an  array  of 
normal-displacement  structural  shape  functions.  The  result  can 
be  expressed  as 

bs  =  pLsxs,  (2.57) 

where  the  matrix  Ls  is  given  by 

Ls=  NNg  Ts  dB,  (2.58) 

JBS 

in  which  Ts  is  a  diagonal  matrix  of  normal  direction  cosines. 

In  practice,  the  entries  of  L8  need  not  be  explicitly  calculated 
and  stored.  Instead,  the  whole  process  of  going  from  vector  x® 
to  vector  bs  can  be  conveniently  packaged  within  a  numerical 
integration  framework.  The  effective  result  of  the  numerical  in¬ 
tegration  process  can  be  presented  in  the  "lumped  area"  form 

bs  =  pAsG‘sxs,  (2.59) 

where  As  and  Gp  are  the  same  as  in  Equation  (2.53). 

Staggered  Integration.  The  semi-discrete  equations  of  motion  of 
the  two  interacting  fields:  submerged  structure  and  fluid  volume, 
are  numerically  integrated  with  a  staggered  solution  procedure  in 
which  only  two  vectors:  nodal  pressures  p  and  structure  displace¬ 
ments  xs,  are  passed  back  and  forth  between  the  structure  and 
fluid-volume  software  modules. 

The  structural  equations  of  motion  (2.52)  are  treated  by  an 
implicit  time  integration  formula,  which  yields  a  (generally  non¬ 
linear)  algebraic  system  that  must  be  solved  at  each  step: 


Esx*  =  h*  +  f?  +  6%, 


(2.60) 
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where 


Es  =  Ms  +  SCs  +  S2Ks,  (2.61) 

6  is  an  integrator-dependent  generalized  stepsize,  h®  is  a  term  that 
embodies  the  effect  of  previous  structural  solutions,  and  t„  is  a 
nonlinear  pseudo-force  term. 

Combining  the  previous  equations  with  the  artificially  damped 
central-difference  advancing  step  (2.35)-(2.39),  the  following  time 
marching  scheme  results: 


Esxs„  =  h*  +f~— «2GsAsPn,  (2.62) 

bs„  =  pAsG‘<,  (2-63) 

§n  =  (Bn  ~  Sn-i)/h,  (2.64) 

$n+ 1/2  =  1/2  +  h  (p„  -  p^  +  (3hc2  Sn),  (2.65) 

^n  + 1  —  ^ n  "T  ^  ^n  +  i/2)  (2.66) 

Q  sn+1  =  —  H  $n+1  +  b®n+1  +  bdn+l,  (2.67) 

Pn  +  i  =  max(pw  +  c2sn+1,0).  (2.68) 


The  only  undefined  term  in  these  equations  is  now  bd,  which 
comes  from  the  DAA  boundary  interaction.  This  term  is  dealt 
with  in  §2.6. 

The  simplest  staggered  solution  procedure  for  the  preceding 
equations  is  obtained  if  one  identifies  b®  +  1  with  b® ,  and  similarly 
for  the  solution-dependent  portion  of  bd.  The  net  effect  of  this 
“last  solution”  staggering  is  that  the  structure  “lags”  one  step 
behind  the  fluid  volume. 

An  obvious  refinement  to  the  previous  scheme  is  the  use  of 
a  predictor  x®p  for  x®  +1  in  (2.63)  instead  of  simply  inserting  the 
last  solution  x® .  For  example, 

K  =  pAs G‘  xf  =  PA, G<  (x'  +  «*).  (2.69) 

Predictors  may  be  used  not  only  to  improve  accuracy,  but  also 
stability  of  specific  integration  formulas ,  as  shown  in  the  following 
study. 
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Stability  of  Staggered  Integration.  The  stability  analysis  presented  here  fol¬ 
lows  Fourier  techniques  similar  to  those  used  in  §2.4.  Let  x  be  the  amplitude 
of  a  structural  displacement  mode,  and  let  m,  a  and  /  be  the  corresponding 
generalized  mass,  contact  area,  and  pressure  force,  respectively.  The  struc¬ 
tural  damping  can  be  neglected  from  the  outset.  The  structure  stiffness 
can  also  be  ignored,  as  a  deeper  analysis  (not  reported  here)  shows.  Stiffness 
nonlinearities  can  be  therefore  ignored.  Thus  only  mass  and  contact  area 
govern  stability. 

The  modal  structural  equations  will  be  integrated  by  the  general  one- 
step  implicit  method: 


Xn  =  Ifl-1  +  h.[tpxn  +  (1  —  <p)xn-i],  (2.70) 

in  =  in-1  +  b[<pxn  4-  (1  —  ^)i„_i],  (2.71) 

This  integrator  specializes  to  Backward  Euler  for  (p  =  1  and  to  the 
trapezoidal  rule  for  <p  =  1/2.  The  generalized  stepsize  is  6  =  tph.  For  zero 
damping  and  stiffness,  the  implicit  equation  (2.60)  reduces  to 

mxn  =  m(xn- 1  -j-  hXn- i)  +  h?[tp2fn  +  ip{  1  —  tp)fn- 1|.  (2.72) 

The  assumed  predictor  for  the  interaction  term  is 

x„+1  =  Xn  +  a  hxn,  (2.73) 

where  a  is  a  free  parameter. 

As  for  the  fluid-volume  modal  equations,  (2.41)  through  (2.44)  apply  with 
only  the  following  changes:  7  is  set  to  one  (pressurized  fluid,  as  a  cavitating 
fluid  mode  maintains  zero  pressure  and  does  not  interact),  and  the  structure 
boundary  coupling  term  bs  is  added  to  (2.43).  This  term  is  divided  by  q, 
which  is  the  generalized  capacitance  associated  with  the  fluid- volume  mode 
under  consideration.  Here  is  the  complete  set  of  difference  equations: 


m(xn_  j 

+  bxn 

-1)  —  ah2\<p2pn  4-  <p(  1  —  <p)pn-i\, 

(2.74) 

K 

=  pa{xn  4-  ahin), 

(2.75) 

V'nH M/2 

0 

=  'fin- 

-1/2  +  MP"  4-  0C2{sn  ~  *n  — 1)]> 

(2.76) 

V'n-f 

•1  =  ^  +hrpn  +  1/2, 

(2.77) 

=  —\ipn+i  4-  K/<I . 

(2.78) 

Pn  + 1  =  C^Sn+1- 

(2.79) 
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Elimination  of  the  intermediate  fluid  variables  s  and  p  furnishes  two 
coupled  difference  equations: 

V'n+i  —  (2  —  £(1  +  P)}lpn  +  (1  —  tPHn-l  = 

{pah,2c2/q)[{l  4-  p)xn- 1  —  pxn-2  +  a(l  +  P)hxn-i  —  <xphxn-2),  (2.80) 
Xn  —  (1  —  <p2fi)xn-  1  +<p(  1  —  tp)n  Xn-2  —  (1  —  Ct(i)hXn- 1 
+^(1  —  ip)apLh2hxn-2  =  $a/m[<p2rpn  +  —  (p)ipn- i],  (2.81) 

where  f  is  given  by  (2.46),  and 


pa2c2h 2 

M  =  - 

mq 


(2.82) 


is  a  dimensionless  parameter  that  measures  the  strength  of  fluid-structure 
modal  coupling. 

Direct  elimination  of  the  velocity  terms  in  the  preceding  equations  is 
messy.  It  is  more  convenient  to  pass  to  the  transform  space  z  first,  and  then 
eliminate  them  through  the  operator  relation 


h°Xk 


- - - Xk,  k  —  n,n  —  1 , . . . 

ipz  4*  1  —  <P 


(2.83) 


*  :  i  *  i »  *  a  *■ 

Oiiiiil  ivmiww.w  Hum  M  n  iidi  v>i  in  i  j|^ 

nomial  can  be  expressed  as 


f  r\  >—  r\\ 

1*.  i 


f' \  ,  - - u: - ^ 

i  lie  icouti/iiig  UKtidticiiovlt 


C{z)  =  Cf/(z)Css(z)  -  KCgj{z)CJe{z),  (2.84) 

where 

CgK(z)  =  23  4-  (—2  4-  {<p2  4-  a <p)pi\z2 
4-  {(1  4-  [2<p(l  —  <p)  4-  <*(1  —  2<p)]pL  }  *4-  (1  —  <p)(l  —  <p  —  a)n,  (2.85) 
Cff(z)  =  {z2~[ 2  -  ?(1  4-  P))z  (2.86) 

C/8(z)  =  (V>4-«)(H-^)2r2-h[(l-|-/&— (V?4-a)(H-2^)]2:— /&(!— a),  (2.87) 

CSf(z)  =  (ipz  4-  1  —  ip).  (2.88) 


Observe  that  C//(z)  and  Css(2)  with  \i  =  0  are  the  characteristic  poly¬ 
nomials  for  the  uncoupled  fluid-volume  and  structure,  respectively,  while 
Cfs{z )  and  C8f(z)  account  for  the  cross-coupling.  Plainly  the  stability  region 
of  C(z)  cannot  extent  beyond  that  of  the  uncoupled  components.  Since  CfS(z) 
is  absolutely  stable  (A-stable)  with  the  choice  1/2  <  ip  <  1,  the  stability 
region  of  (2.84)  must  lie  inside  that,  of  Cj j(z),  which  is  in  fact  given  by  (2.50). 

A  long  but  straightforward  analysis  shows  that  the  quintic  polynomial 
C(z)  has  a  double  root  at  2=1.  Removing  the  factor  (z —  l)2  reduces  C(z) 
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to  a  cubic  polynomial: 


C(z)  =  z  Cj(z )  —  OJLCm(z), 
where  £  is  given  by  (2.46),  and 

Cf(z)  =  a2  -  [2  -  f(l  +  0)]z  + 


(2.89) 


(2.90) 


Cm(z )  =  <p(<p  +  a)  z2  +  2<p(l  —<p)  +  a(l  —  2<p)z  +  (1  —  v>)(l  —  <P  —  a) 
=  (<pz  +  1  —  <p)\((p  +  a)z  +  1  —  <p  —  a)],  (2.91) 


[Note  that  Cr(z)  is  precisely  (2.47).)  The  stability  of  C(z)  was  studied  with  a 
computer  program.  Some  results  of  the  study  are  shown  in  Figures  2  through 
5.  In  these  figures  the  stability  region  is  plotted  in  the  h,0  plane  over  the 
"window" 

0  <  h/hc  <1.  0  <  0  <  1,  (2.92) 


where  he  is  the  critical  (Courant)  timestep  (2.32)  for  an  uncoupled  and  un¬ 
damped  (0  =  0)  fluid-volume  mesh.  The  other  three  free  parameters  are  <p, 
a  and  m.  which  characterize  the  integration  formula  (2.70)-(2.71),  the  inter¬ 
action  term  predictor  (2.75).  and  the  modal  coupling  strength,  respectively. 
The  parameter  f  is  implicitly  defined  by  h/hc .  since  f  =  A{h/hc )2  from 


(2.19)  for  0  —  0  and  Tf  =  1 . 

Figure  2  show  stability  regions  for  the  trapezoidal  rule  (<p  =  1/2)  when 
a  =  0,  i.e.  the  last  displacement  solution  is  used  in  the  interaction  term  6s. 
as  in  (2.t>3).  Tiie  eight  fr/mc  show  n  in  Figure  2  pertain  to  fixed  values  of 


Mr  = 


pa2c'2h ^ 
mq 


(2.93) 


which  is  simply  (2.82)  evaluated  at  the  Courant  timestep.  This  is  the  critical 
physical  parameter  as  regards  stability  of  the  staggered  time-integration 
procedure.  The  values  of  Mr  for  the  eight  frames  are  listed  in  the  figure 
caption.  Stable  regions  are  dark  shaded.  The  stable  region  for  mc  =  0 
corresponds  to  the  no-interaction  case,  whose  equation  is  (2.50).  It  can  be 
seen  that  as  mc  increases,  the  fluid  damping  coefficient  p  can  have  a  dramatic 
eflcct  on  stability.  For  example,  if  mc  =  1>  the  largest  stable  h  is  virtually 
zero  if  0  =  0,  but  surges  to  about  0.81  he  if  0  =  0.25.  This  effect  should 
be  contrasted  to  the  uncoupled  case  studied  in  §2.4,  in  which  increasing  0 
always  reduces  the  stable  stepsize.  But  even  with  damping  help  in  the  range 
0  <  0  <  1.  no  stable  stepsize  in  the  window  (2.91)  essentially  survives  for 
Me  >  5. 

Figure  3  also  pertains  to  the  trapezoidal  rule,  but  now  a  is  1,  which 
effectively  amounts  to  using  the  predictor  (2.69).  It  can  be  seen  that  use  of 
this  predictor  substantially  extends  the  stability  region  for  large  Me-  The 
effect  of  the  damping  coefficient  0  is  not  so  dramatic  as  in  the  previous  case, 
but  it  is  still  pronounced  for  high  Mr  values. 
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“Overpredicting"  with  a  >  1  has  even  more  beneficial  effects  on  stability, 
as  evidenced  by  Figure  4,  which  corresponds  to  a  =  2.  But  accuracy  suffers, 
so  using  a  =  1  is  recommended. 

Finally,  Figure  5  shows  results  for  the  Backward  Euler  method  {<p  =  1). 
This  integrator  introduces  high  numerical  damping  in  the  structure,  whereas 
the  trapezoidal  rule  introduces  none.  The  stability  regions  for  large  hc  are 
substantially  enlarged,  and  are  now  fairly  insensitive  to  both  {}  and  a  (effect 
of  the  latter  is  not  shown  here.) 

Mesh  Sizing  Guidelines.  The  fact  that  hc  is  the  critical  stability  parameter 
for  fluid-structure  coupling  may  be  used  to  derive  some  dimensioning  guide¬ 
lines  for  the  fluid-volume  mesh. 

It  is  assumed  that  the  wet-surface  structure  is  a  shell,  the  discretization 
of  which  is  known  a  priori.  Now  consider  the  interaction  between  two  ad¬ 
jacent  physical  elements:  (1)  a  square  plate  dimensioned  L  X  L,  with  thick¬ 
ness  ts  and  density  pg,  and  (2)  a  rectangular  fluid-volume  brick  dimensioned 
Lx  LxD,  where  D  <  L  is  the  dimension  normal  to  the  plate.  The  contact 
area  is  L  X  L.  Replacing 

h2c=D2/c2,  a2  =  L4,  m  =  p,L2ts  q  =  L2D/2,  (2.94) 


into  (2.93)  yields 


Me  = 


2  pD 
Pstg 


For  sieei  in  wiier,  p.  /p  8,  so  ihai 


He  PH 


4t, 


(2.95) 


(2.96) 


This  can  be  used  for  an  order-of-magnitude  (generally  conservative)  estimate 
of  nc  high-frequency,  localized  “mesh  modes”.  The  estimate  helps  in 
sizing  the  first  layer  of  fluid-volume  elements  adjacent  to  the  structure.  The 
mesh  can  then  be  “radially"  continued  into  the  fluid  as  far  as  necessary. 

For  example,  suppose  that  ts  =  2  inches,  and  that  it  is  desired  to  keep 
He  <  5  from  stability  considerations.  Then  D  should  not  exceed  40  inches. 

This  simple  rule  has  been  used  to  size  the  fluid-volume  meshes  for  the 
example  problems  presented  in  Section  IV,  so  that  He  does  not  exceed  5. 
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Figure  2.  Stability  of  (2.84)  for  trapezoidal  rule  (<p  —  0.5) 
and  last-solution  predictor  a  —  0.  Each  frame  covers  the 
“window"  0  <  h/hc  <  1  horizontally  and  0  <  (3  <  1 
vertically.  Stable  regions  are  dark  shaded.  Starting  from  the 
upper  left  corner,  the  eight  frames  correspond  to  the  following 
values  of  nc'-  0,  0.1,  0.5,  1.0,  2.0,  5.0,  10.0  and  20.0. 


Figure  3.  Stability  of  (2.84)  for  trapezoidal  rule  (ip  =  0.5)  and 
full-step  predictor  a  =  1.  Each  frame  covers  the  “window” 
0  <  h/ he  <  1  horizontally  and  0  <  f)  <  1  vertically.  Stable 
regions  are  dark  shaded.  Starting  from  the  upper  left  corner, 
the  eight  frames  correspond  to  the  following  values  of  pc-  0, 
0.1,  0.5,  1.0,  2.0,  5.0,  10.0  and  20.0. 


Figure  4.  Stability  of  (2.84)  for  trapezoidal  rule  (<p  =  0.5) 
and  “overextrapolated”  predictor  a  =  2.  Each  frame  covers 
the  “window”  0  <  h/hc  <  1  horizontally  and  0  <  p  <  1 
vertically.  Stable  regions  are  dark  shaded.  Starting  from  the 
upper  left  corner,  the  eight  frames  correspond  to  the  following 
values  of  /ic:  0,  0.1,  0.5,  1.0,  2.0,  5.0,  10.0  and  20.0. 
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Figure  5.  Stability  of  (2.84)  for  Backward  Euler  (<p  =  1) 
and  last-solution  predictor  a  =  0.  Each  frame  covers  the 
window"  u  <  H/hc  <  l  horizontally  and  0  <  p  1 
vertically.  Stable  regions  are  dark  shaded.  Starting  from  the 
upper  left  corner,  the  eight  frames  correspond  to  the  following 
values  of  uc:  0,  0.1,  0.5,  1.0,  2.0,  5.0,  10.0  and  20.0. 


§2.6  Fluid- DAA  Interaction 

Role  of  DAA  Boundary.  The  Doubly  Asymptotic  Approximation 
(DAA)  boundary  truncates  the  fluid-volume  mesh  to  finite  extent. 
In  its  discrete  form,  it  consists  of  boundary  elements  in  contact 
with  faces  of  fluid-volume  brick  elements.  This  computational 
field  should  ideally  operate  as  a  transparent  entry  boundary  for 
incoming  (incident)  waves,  and  as  a  perfectly  radiating  boundary 
for  outgoing  waves.  Because  of  the  nature  of  the  DAA,  these  con¬ 
ditions  are  asymptotically  satisfied  in  the  limit  of  high-frequency 
and  low-frequency  motions,  and  approximately  otherwise. 

The  DAA  boundary  element  mesh  is  usually  constructed  so 
that  its  nodes  coincide  with  fluid- volume  nodes.  The  net  result 
is  that  DAA  elements  lie  on  brick  faces.  But  all  DAA  computa¬ 
tional  vectors  are  expressed  not  in  terms  of  nodal  point  values, 
but  rather  of  values  at  control  points,  which  are  located  at  the 
centroid  of  each  DAA  element.  So  it  sometimes  becomes  neces¬ 
sary  to  distinguish  quantities  such  as  displacement  and  pressure 
vectors,  area  matrices,  etc.,  which  can  be  referred  to  either  set  of 
points. 

In  this  and  following  sections,  letters  /  and  d  applied  as 
subscripts  or  supercripts  to  a  matrix  or  vector  symbol  are  used  to 
indicate  that  it  pertains  to  DAA  cont  rol  points  and  to  fluid-volume 
nodes  located  on  the  DAA  boundary,  respectively;  for  example, 
total-pressure  vectors  and  pd.  If  neither  appears,  d  is  assumed. 

Boundary  Interaction  Terms.  The  DAA  boundary  acts  on  the  fluid 
volume  through  the  forcing  term  bd  of  (2.24).  Arguments  similar 
to  those  used  in  §2.5  can  be  offered  to  derive  a  formal  matrix 
expression  that  relates  brf  to  the  vector  of  global  (X,Y,Z) 
displacements  at  the  DAA  control  points: 

(2-97) 

Here  A d  is  a  diagonal  matrix  of  contributing  areas  that  surround 
fluid-volume  nodes,  Gdf  is  a  transformation  matrix  from  DAA 
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control  points  to  fluid- volume  nodes,  and  T/  is  a  diagonal  matrix 
of  direction  cosines  of  the  boundary  normal  (positive  going  into 
the  fluid)  evaluated  at  the  DAA  control  points.  As  in  the  case 
of  the  structure  interaction,  entries  of  A d  and  Gdf  need  never  be 
explicitly  computed,  for  the  whole  vector  transformation  process 
is  elegantly  hidden  by  one-point  Gauss  isoparametric  integration. 

The  DAA  displacement  vector  may  be  decomposed  into 
three  components  due  to  the  free-field  incident  wave,  scattered 
wave,  and  hydrostatic  pressure: 

x/  =  xfi  _|_  xfs  _|_  xfH  (2.98) 

The  hydrostatic  displacement  (but  not  the  pressure)  may  be 
set  to  zero  ab  initio ,  as  it  cancels  out  in  the  relative  displacement 
formulation  used  here.  The  other  two  components  are  studied  in 
the  following  subsections. 


Incident  Wave.  We  consider  incident  spherical  and  plane  waves. 
An  incident  snherical  waveform  is  comDletelv  defined  bv  eivine 
its  origin  (charge  location)  and  the  pressure  profile 

p'(t  —  t0:R)  (2.99) 

recorded  at  a  reference  location  whose  distance  to  the  wave  origin 
is  R.  As  the  wave  clock  can  be  adjusted  through  an  arbitrary  time 
shift  t0 ,  the  reference  location  may  be  conveniently  specified  as 
the  fluid-volume  node  “touched”  by  the  w  avefront  at  the  reference 
time  t  —  0;  for  this  node  R  —  Ro.  The  free-field  incident  pressure 
can  then  be  readily  calculated  at  any  fluid  point  for  all  times 
—  Rq/c  jC  t  <C  -}-oc. 

The  free-field  fluid-particle  displacement  at  an  arbitrary 
point  joined  to  the  w'ave  origin  through  the  u.iit  direction  vector 
ft  may  be  determined  from  the  relation 

y'u)  =  />*')/!  (2.ioo) 

pc  pR 
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where  each  superscript  asterisk  denotes  temporal  integration  from 
t  =  —Ro/c  through  t.  This  formula  can  be  specialized  to  DAA 
control  points  to  provide  the  matrix  expression 

T/x/7  =  —  +  -R_1r/  V/7,  (2.101) 

pc  p 

where  R  is  a  diagonal  matrix  containing  distances  from  the  wave 
origin  to  DAA  control  points,  and  Tf  is  a  diagonal  matrix  of  the 
cosines  of  the  angles  between  the  local  propagation  direction  ft 
and  the  outward  unit  normal  vectors  at  the  DAA  control  points. 

For  an  incident  plane- wave,  R  *-*■  oo,  and  the  second  term  on 
the  right  of  (2.101)  drops  out. 

Scattered  Waves.  Displacements  caused  by  scattered  waves  are 
calculated  from  the  simplest  Doubly  Asymptotic  Approximation, 
which  reads 


MfpfS-j-  pcAfpfS  =  pcMf(ifS.  (2.102) 

In  this  equation,  p/s  and  u^s  are  column  vectors  of  scattered- 
wave  pressures  and  normal  fluid-particle  velocities,  respectively, 
at  DAA  control  points,  M/  is  the  (fully  populated)  mass  matrix 
for  irrotational  incompressible  motions  of  the  fluid  external  to  the 
DAA,  and  A/  is  a  diagonal  matrix  of  boundary  element  areas. 

To  get  x-fb,  integrate  (2.102)  twice  in  time,  and  solve  for  the 
scattered  normal  displacements  that  appear  on  the  right-hand 
side: 

T/  xfs  =  (2.103) 

where  constants  of  integration  for  and  are  determined 
from  the  initial  conditions  discussed  in  §2.7. 

Remark.  Do  not  confuse  A/  with  the  A d  that  appears  in  (2.97).  Both  are 
diagonal  area  matrices  but  the  first  one  pertains  to  DAA  control  points  and 
the  second  one  to  fluid-volume  nodes  at  the  DAA  boundary. 


2-29 


Mass/Dampiag  Split.  The  two  components  of  x^s  that  appear 
in  (2.103)  have  different  physical  significance  and  deserve  to  be 
identified  separately: 


xfD  =  irr1  f>fs,  (2.104) 

pc  1 

x/M  =  T Jl  My1A/  (2.105) 

Here  superscripts  D  and  M  stand  for  damping  and  mass,  respec¬ 
tively,  in  accordance  with  the  following  interpretation. 

The  displacement  vector  x^D  corresponds  to  the  DAA  operat¬ 
ing  as  a  “pc  boundary”,  radiating  high-frequency  energy  out  into 
the  external  fluid.  The  displacement  vector  x^M  corresponds 
to  the  DAA  acting  as  an  “added-mass  boundary”,  accounting 
rigorously  for  low-frequency  “sloshing”  of  the  external  fluid.  As 
discussed  later,  these  different  interpretations  translate  into  dif¬ 
ferent  numerical  treatments  in  the  implementation  of  the  time 
integration  procedure. 

Insertion  of  the  various  displacement  terms  into  (2.97)  — 
with  the  hydrostatic  component  excluded  —  splits  the  interaction 
vector  into  three  components: 


bd  =  b1  +  bD  +  bM 

(2.106) 

b^ArfG^fV  +  R-'^'), 

(2.107) 

bD  =  -AdGdff>-rs, 

C 

(2.108) 

bM  =  pAdGdfMj'Af)t's. 

(2.109) 
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Scattered  Pressure  Calculations.  The  fluid  volume  acts  on  DAA 
boundary  elements  through  the  total-pressure  vector  p,  which  is 
obtained  in  the  course  of  the  time  integration  solution  process. 
Total  pressures  mav  be  interpolated  from  fluid-volume  nodes  to 
DAA  control  points  through  the  transformation  matrix  G/d  = 


G  W- 


y'  =  G 


(2.110) 


The  scattered  pressure  component  now  follows  by  substracting  off 
the  incident  and  hydrostatic  components: 


p/5  _  p /  _  p //  _  p /// 


(2.111) 


This  scattered  pressure  data  can  be  time-integrated  numeri¬ 
cally  to  generate  jj yf^  and  *pfs.  For  example,  using  the  trape¬ 
zoidal  rule: 


H'l,  =  H*  +  2{pfnS  + 


(2.112) 


(2.113) 


Finally,  these  pressure-integral  values  can  be  inserted  into 
(2.108)  and  (2.109)  to  close  the  interaction  loop. 


Pressure  Correction.  Computational  experiments  with  a  fully 
staggered  solut  ion  procedure  for  the  DAA  interaction  have  shown 
that  undesirable  pressure  oscillations  develop  near  the  DAA  boun¬ 
dary.  These  oscillations  are  caused  by  the  time  lag  in  the  treat¬ 
ment  of  the  “pc- boundary”  coupling  term  (2.104).  This  lag  in¬ 
terferes  with  the  energy-radiation  process  for  outgoing  scattered 
waves.  The  spurious  pressure  oscillations  eventually  reflect  back 
to  the  structure  and  distort  its  response. 

The  problem  has  been  solved  by  using  a  simultaneous  pres¬ 
sure  solution  for  the  /?r-boundary  term,  while  the  “added  mass” 
term  (2.105)  is  treated  by  staggered-solution  techniques.  The 
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whole  business  can  be  transacted  at  the  node  level.  Thus,  consider 
an  individual  fluid-volume  node  located  at  the  DAA  boundary. 
Rewrite  the  interaction  term  (2.108)  as 

bD  =  lad/c)}s,  (2.114) 

where  bD ,  ad  and  f)s  denote  entries  of  vectors  bD,  diag(Ad)  and 
G  respectively,  pertaining  to  the  node  under  consideration. 

Next,  time-discretize  (2.114)  through  the  trapezoidal  rule  (2.112): 

bn  +  !  +{h/2  ){Pn+Pn  +  i)}-  (2.115) 

Now  use  p*+l  =  p„+l  -  p[n+l  -  pjf+1  to  get 

bn  + 1  =  (ad/c)lf>Sn  -\-{h/2)(p%  -  pln+l  -  Pn  +  1+Pn  + 1 )] -  (2-116) 

According  to  the  problem-modelling  assumptions  stated  in 
§2.1,  cavitation  should  not  occur  at  the  DAA  boundary.  Hence 
the  nodal  pressure  calculation,  given  by  (2.67)-(2.68),  becomes 

7 (Pn+.  ~p")  =  c2(  r  +  />'„+,  +  b™)  + 

+  ^(P'l  -  pi+i  -  p"+,  +  P»+.)]  (2.117) 

where  r,  q,  b1  and  6iU  are  entries  of  vectors  r  =  —  H^,  diag(Q),  b; 
and  bA/,  respectively,  for  the  node  under  consideration.  The  next 
total-pressure  value  pn+\  appears  on  both  sides  of  this  equation. 
For  simultaneous  solution,  the  unknown  term  is  moved  to  the  left, 
giving 

qPn+i  =  qil+K)pn+i  =  g,  (2.118) 

where 

k  =  —chad,  (2.119) 

and  g  embodies  all  “leftover”  right-hand  side  components.  The 
net  effect,  of  all  this  is  that  diagonal  entries  of  Q  pertaining  to 
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fluid-volume  nodes  on  the  DAA  boundary  must  be  multiplied  by 
a  correction  factor  1  -f-  k. 

Remark  l.  The  correction  factor  (1  +  /c)  is  stepsize  dependent,  and  tends  to 
one  as  h  *-*•  0.  It  is  fairly  easy  to  show  that  k  <  1  if  h  <  he,  and  takes  the 
value  1  on  a  regular  grid  if  h  equals  the  Courant  stepsize  he- 

Remark  2.  Kq  is  half  of  the  fluid  volume  swept  by  an  area  a*  over  the  distance 
ch  travelled  by  a  sound  wave  over  the  time  increment  h. 

Remark  3.  Since  the  pressure  correction  is  node-level,  the  basic  philosophy 
of  the  staggered  solution  procedure,  which  calls  for  only  vector  transfer 
information,  is  not  violated.  Had  the  simultaneous  solution  procedure  been 
extended  to  include  the  fluid-mass  interaction  vector  bM,  the  fluid  volume 
analyzer  would  need  to  know  about  the  full  matrix  M /. 

Flemark  4.  In  the  terminology  of  coupled-system  partitioned  analysis,  the 
process  by  which  selected  field  quantities  are  manipulated  into  the  left-hand 
side  of  the  equations  of  another  field  is  called  augmentation.  This  technique 
is  primarily  used  to  improve  numerical  stability  characteristics  [8]. 

Staged  DAA  Analysis.  Introduction  of  the  pressure  correction 
mechanism  effectively  splits  the  analysis  of  the  DAA  field  into 
two  stages.  In  the  first  stage,  the  DA\  displacement  vector 

x/  =  *{+i  +  *lD  +  xi+f  (2.120) 

is  evaluated  and  supplied  to  the  fluid  volume  analyzer.  In  (2.120), 
+  J  is  a  predictor  for  x{^j.  This  term  is  generated  by  ex¬ 
trapolating  the  double  integral  of  scattered  pressure,  for  example 

V&?  =  ViS  +  Af£S,  (2-121) 

which  is  then  inserted  in  (2.109).  The  primary  advantage  of 
predicting  b{+,,  rather  than  using  the  previous  value,  is  to  im¬ 
prove  numerical  stability  characteristics,  as  shown  later  in  this 
section. 

The  second  stage  begins  on  exit  from  the  fluid-volume  ana¬ 
lyzer,  which  returns  the  total  pressure  vector  pj|  .  Now  the 
scattered  pressure-integral  terras  can  be  corrected  using  Equations 
(2.110)  to  (2.113). 
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For  the  implementation  of  this  staged  process  it  is  more 
convenient  to  use  an  incremental  formulation  based  upon  the  total 
pressure  rather  than  the  scattered  pressure,  as  described  next. 

Incremental  Formulation.  Although  the  value  of  accumulated 
since  t  =  0  given  in  (2.120)  must  be  supplied  to  the  fluid- volume 
analyzer,  in  practice  it  is  convenient  to  handle  the  DAA  computa¬ 
tions  incrementally ,  since  some  terms  are  known  precisely  while 
others  must  be  estimated  and  then  subsequently  corrected. 

The  known  terms  involve  the  incident  and  hydrostatic  pres¬ 
sure.  Some  of  them  arise  because  of  the  DAA  formulation,  which 
is  based  on  the  scattered  pressure,  and  the  need  to  work  with 
total  pressure  for  the  implicit  treatment  of  the  pc  boundary  term. 
They  are: 

A  *"■«+,-«'.  (2.122) 

A?"  =  (fc/2XK'+1  +  H'),  (2-123) 

Ap/W  =  hpH ,  (2.124) 

K2.  =  K"  +  A^H,  (2.125) 

AVVH  =  A(l  +  A)J/ffi.  (2.126) 

Hence 

T /  A xinown  =  -T(A i>f'  +  A^H)  -  Mj'A/fAy"  +  aY/H) 

+  (2.127) 

pc  p 

The  unknown  terms  that  must  be  estimated  involve  integrals  of 
the  total  pressure  as 

aF  =  hp{,  (2.128) 

ipf  =  l>hl,  (2.129) 
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(2.130) 


T;  Ax{s„m„(,  „  =  1-Ay  +  M-f'Af  Ap*f- 

The  predicted  value  of  x/t_M  that  is  supplied  to  the  fluid- volume 
aualyzer  is  then 

^*{:nown  estimated •  (2.131) 

On  return  from  the  fluid-volume  analysis,  displacements  and  pres¬ 
sures  are  corrected: 


aF  — </>/2Hp{+i  +  p£)> 

(2.132) 

H  +  ,  =  Pn  +  Ap7, 

(2.133) 

AV/  =  (/!/2)(K+,+^), 

(2.134) 

(2.135) 

known  corrected ’ 

(2.136) 

Notice  that  in  this  formulation  there  is  no  need  to  keep  track 
of  the  accumulated  double  integral  of  pressure 

Stability  of  Staggered  Integration.  The  following  stability  analysis  only  in¬ 
vestigates  the  effect  of  staggering  the  fluid-mass  coupling  term  bM.  The 
radiation-damping  term  bD  has  no  effect  on  stability,  for  it  is  treated  im¬ 
plicitly.  Incident  and  hydrostatic  components  are  dropped  and  scattered 
pressure  becomes  pressure.  Cavitation  is  ignored.  Absolute,  rather  than  in¬ 
cremental,  quantities  are  used.  Under  these  assumptions,  the  modal  equations 


for  t he  advancing  slop  read: 

$n+  1/2  =  1/2  +  Ml1  +  P)Pn  ~  PPn- l],  (2.137) 

tpn  +  l  =  Ifin  +  (2.138) 

V£+1=V.  +«*»-.  (2.139) 

b*1  =P«dt*P*Z+ 1.  12.140) 

Sn  +  i  =  -f  b*f/q,  (2.141) 
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2 

pn-t-1  —  C  Sri+1* 


$n  +  l  =  "E  (V^)  (Pn  "H  Pn  +  l)> 
*P*„-f  i  =  *P*n  +  (V2)  (K  +  £n  +  l)- 


(2.142) 

(2.143) 

(2.144) 


Io  these  equations,  ip,  s  and  p  are  amplitudes  of  'l',  s  and  p,  respectively,  for  a 
fluid-volume  normal  motion  of  eigenvalue  X  as  in  §2.5;  ad  and  q  are  general¬ 
ized  values  of  A d  and  Q,  respectively,  the  latter  embodying  appropriate  pres¬ 
sure  correction  factors  (2.119);  O!  is  a  pressure-integral  predictor  coefficient; 
and  (  are  the  roots  of  the  “fluid  boundary  mode”  symmetric  eigenproblem 


£  M/  w  =  A/  w, 


(2.145) 


w  being  the  boundary  mode  excited  by  the  modal  volume  pressure. 

Elimination  of  the  intermediate  variables  ip  and  s  yields  two  coupled 
difference  equations: 


1pn+X  —  2 1pn  +  =  — fl2[(l  +  P)pn  ~  Ppn- i] 

P»  +  i  —  4ish~2(*p*n  -f  OLh^n)  =  —  Xc2^n  +  1, 


in  which 


pc2ad  ( 


(2.146) 

(2.147) 


(2.148) 


is  a  dimensionless  modal-couDling  coefficient  that  Dlavs  a  role  analogous  to 
that  of  fj,  in  §2.5.  (The  factor  4  is  introduced  for  convenience  in  subsequent 
manipulations.)  Transform  these  difference  equations  to  the  z  plane,  and 
eliminate  the  pressure-integral  terms  through  the  operator  relations 


**  __  h  z  -f  1  *  _  h2 ( 2+ IV 
V  n  —  ~  7  Pn  —  —1  “  1  Pn 

2  z  —  1  4  \z  —  1  / 


The  resulting  cubic  characteristic  polynomial  is 

C(z )  =  zCj(z)  —  ci 'cCm(z), 
where  f  is  given  by  (2.46), 

Cf(z )  =  z2  —  (  2  —  C(1  +  P)\z  -hi  —  f/9, 

Crrriz)  =  {Z  +  l)2  +  2C*(*2  —  1), 
and  vc  is  v  evaluated  at  the  Courant  timestep: 


(2.149) 


(2.150) 


(2.151) 

(2.152) 


Vc  = 


pc2adZhl 


(2.153) 
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The  similarity  of  these  expressions  with  (2.89)  through  (2.91)  is  apparent. 
Polynomial  C/(z)  is  the  same  as  (2.90),  and  Cm(a)  is  precisely  (2.91)  if  <p  = 
1/2.  Thus  the  stability  regions  of  (2.89)  for  the  trapezoidal  rule  and  of  (2.151) 
coincide  if  pic  and  vc  are  identified.  One  can  therefore  refer  to  Figures  2 
through  1  for  the  stability  of  C(z )  with  varying  a,  simply  by  replacing  ptc 
by  vc . 


Estimating  vc.  To  apply  the  conclusions  of  the  stability  analysis,  it  remains 
to  obtain  a  ballpark  upper-bound  estimate  for  vc ■  To  get  it,  an  admittedly 
idealized  situation  is  considered.  Imagine  that  the  DAA  boundary  is  a  sphere 
of  radius  R.  and  that  adjacent  fluid-volume  elements  are  (roughly)  bricks 
dimensioned  L  X  L  X  L.  Assume  further  that  the  fluid  motions  are  axisym- 
metric  fluid- boundary  modes  of  a  sphere  (Legendre  functions)  with  circum¬ 
ferential  wave  number  m.  The  corresponding  eigenroot  (  is  then  given  by 


_  2m-f-l  _  TO  -)-  1 

TO/-  4  k  pR 3  pft 

(m  -f- 1  )(2w  +1 ) 


(2.154) 


Next,  insert  ad  =  L2,  q  =  L3/ 2(1  -(-  ac),  he  =  L/c,  and  (2.154)  into  (2.153) 
to  get 

(I  +K)(m-f  I)L 

^  =  ^ - ■  (2.155) 

As  the  coup  de  grace,  claim  that  the  boundary  mode  wavelength  is  of  the 
order  of  the  fluid-volume  mesh  size  L,  so  that 


r  27 rR 

L  Fzi  - . 

TO  -)-  1 


(2.156) 


this  being  clearly  a  worst-case  scenario.  Inserting  (2.156)  in  (2.155)  yields 
Vi -  fx  7r/(  1  -f-  ac).  Since  0  <  k,  <  1, 


vc  <  n  fh  3, 


(2.157) 


which  is  the  estimate  sought. 
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§2.7  Response  Calculation  Details 

This  section  meshes  together  the  theoretical  developments  out¬ 
lined  in  §2.3  through  §2.6  with  a  starting  procedure  to  put  together 
a  practical  response  calculation  scheme. 

The  Reference  State.  Starting  a  staggered  solution  procedure 
that  involves  three  computational  fields  and  combines  explicit  and 
implicit  integration  methods  is  a  delicate  task.  If  the  integration 
does  not  start  right,  it  cannot  be  corrected.  Furthermore,  the 
starting  procedure  should  be  independent  of  whether  the  structure 
is  linear  or  nonlinear,  and  be  readily  extendible  to  internal  fluid 
problems,  in  which  the  motion  of  the  structure  boundary  provides 
the  input  excitation. 

To  meet  these  goals,  the  state  at  t  =  0  is  used  as  reference 
state.  Successive  integration  steps  determine  the  deviation  from 
the  reference  state,  rather  than  the  total  state.  For  example, 
the  total  pressure  vector  at  fluid  volume  nodes  is  actually  broken 
down  into  four  components: 

P  =  VH  +  pi  +  (p7  -  Po)  +  Ps  (2.158) 

o 

and  similarly  for  ty,  'k,  s,  etc.  The  pressure  determined  by  the 
central  difference  scheme  is  the  reference  state  deviator 

P7-PJ  +  PS,  (2.159) 

which  added  to  the  reference  pressure 

Po  =  P^  +  Po  (2.160) 

furnishes  the  total  pressure  p.  Extending  this  idea  to  the  full 
coupled  system  leads  to  the  following  solution  procedure. 


Initialization.  To  start  up  the  time  integration  process,  do  the 
following. 

•  Hydrostatic  Pressure.  Calculate  hydrostatic  pressure  pH  at  all  fluid  volume 
nodes,  structure  wet-surface  nodes,  and  DAA  control  points. 

•  Incident  Pressure.  Given  a  spherical  or  plane  wave  pressure  profile,  origin, 
and  fluid  volume  node  touched  by  the  wavefront  at  t  =  0  (the  “wavefront 
node"),  calculate  the  initial  incident  pressure  Pq  at  fluid  volume  nodes  and 
DAA  control  points  (The  wavefront  node  may  be  in  contact  with  the  struc¬ 
ture,  but  the  front  must  not  intercept  the  structure.) 

•  Reference  Structure  Solution.  Calculate  the  static  response  of  the  structure 
to  the  hydrostatic  pressure;  let  Xq  be  the  corresponding  displacement  vector. 
Initialize  historic  vectors. 

•  Reference  DAA  Solution.  Calculate  fluid-particle  displacements  at  DAA 
control  points  due  to  hydrostatic  pressure  and  initial  incident  wave.  Initialize 
pressure  integral  vectors. 

•  Stepsize.  Select  initial  time  increment  h. 

•  Initial  Velocity  Potential.  Calculate  4'_1/2  at  fluid  volume  nodes  by 
ii  -‘grating  the  incident  wave  flux  from  t  =  — oo  to  t  =  — h/2.  (This  may 
be  done  anahtically  fo,  simple  waveforms,  and  numerically  otherwise.) 

•  Initial  Displacement  Potential.  Calculate  at  fluid  volume  nodes  by 
doubly  integrating  the  incident  wave  flux  from  t  =  — oo  to  t  =  0.  (Same 
remark  as  for  the  velocity  potential  calculation.) 

lime  Integration  Process.  For  n  =  0, 1, . . .,  do  the  following. 

•  Structure  Analysis.  Solve  the  structure  equations  for  the  next  displacement 
vector.  In  th  e  case  of  a  linear  structure, 

E*  (x;  -  x*)  =  K~  62 GsARpr, .  (2.161) 

•  DAA  Boundary  Analysis  (1st  Stage).  The  incremental  formulation  described 

in  §2.6  is  used.  At  DAA  control  points,  compute 

P  l  =G/dpt 

aF'-K'+i-K'. 

AV/;  =  (V2)(K  +  1  +Kf). 
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(2.162) 

(2.163) 

(2.164) 


A$/H  =  hpw,  (2.165) 

H  +  i=HH  +  ±VH,  (2.166) 

A^f//  =h{\  +  ^)H+i-  (2-167) 

T /  A xLown  =  -~(A^;  +  AFH)  -  Mj'AfiA V/J  +  Atf'") 

pc 

+  —  rf  A&1  -h  -r/R^A})*".  (2.168) 

pc  p 

T/  A x'slimaUd  =  -p'  +  (2.169) 

pc 

x{p  =  xf„  +  Ax{nown  +  Ax{stimated.  (2.170) 

•  Fluid  Volume  Analysis.  At  fluid  volume  nodes,  compute 

8„  =  (■„  —tn-i)/h,  (2.171) 

^n-M/2  —  'f'n  — 1/2  ”j“  h  (pn  —  pH  4~  /3h,C2&n)-  (2.172) 

—  'K.  H-  .’i'K+i/o,  (2.173) 

Tn  +  l  =  —  H(*n  +  1  —  *o),  (2.174) 

K  =  pA,Gs(x*  —  xj),  (2.175) 

bf,  =  pArfG(f/T/(x{p  —  Xq  ) ,  (2.176) 

(where  in  (2.171)  8n_!  vanishes  if  n  =  0).  The  pressure  calculation  steps 
depend  on  nodal  boundary  conditions,  and  are  best  stated  in  terms  of  node 
quantities.  For  internal  fluid  nodes, 

Sn+l  =  rn+l/Qj  (2.177) 

Pn+i  =  max(pw  +Po  +  c2sn  +  1,0).  (2.178) 

For  nodes  in  contact  with  the  structure, 

sn  +  l  =  (rn  +  1 +*;)/*  (2.179) 

Pn  +  1  =  rnax(pH  +  pj  +  c2s„  +  1,0).  (2.180) 
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For  nodes  in  contact  with  the  DAA  boundary 

k  =  chad/(2q), 


(2.181) 


>  7»  +  1 


Sn  + 


(r„  +  i  +  b*)/~q  —  sn 


(i  +  #0 

Pn+  1  =  VH  +  C2S„  +  i. 


(2.182) 

(2.183) 


For  nodes  at  a  boundary  of  specified  pressure  p  ( e.g .  a  free  surface), 


Sn  +  l  —  0, 

(2.184) 

Pn-f-1  =  P  ■ 

(2.185) 

•  DAA  Boundary  Analysis  (2nd  Stage).  Return  to  the  DAA  analyzer  to 

correct  control  point  values: 

Pn  4- 1  —  ®/rfPn  +  l- 

(2.186) 

Ap'  =  (ft/2)  (p£h-i  +p£), 

(2.187) 

K+,  =K  +  aF, 

(2.188) 

AV'  =(V2)(K+1  +  K). 

(2.189) 

T/ Corrected  =  ~ +  M/lA/  A*Pf, 
pc 

(2.190) 

X«4-l  ~t~  ^x(notpn  Axcorrec(ed‘ 

(2.191) 

•  Advance.  Increment  counter  n  by  one,  time  t  by  h,  and 
structure  subsystem,  [j 

return  to  the 
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Section  HI 


IMPLEMENTATION  AND  USAGE 

§3.1  Implementation  Overview 

A  number  of  modifications  have  been  made  to  the  standard  USA- 
STAGS  code  in  order  to  accomodate  the  CFA.  In  addition,  one 
important  assumption  has  been  made  that  is  implicit  in  the  treat¬ 
ment  of  Section  II  and  it  should  be  clearly  stated  here.  At  this 
time  the  computational  model  does  not  allow  the  DAA  boun¬ 
dary  to  be  coincident  with  the  structure  boundary  at  any  point. 
In  other  words,  there  must  always  be  at  least  one  layer  of  fluid 
volume  elements  between  the  structure  and  the  DAA  boundary, 
even  if  the  problem  under  investigation  involves  internal  fluid  with 
non-cavitating  external  fluid.  Relaxation  of  this  restriction  will 
be  the  subject  of  a  future  study  into  alternate  forms  of  the  inter¬ 
action  equations. 

The  stability  analyses  and  starting  procedure  described  in 
Section  II  clearly  provide  the  groundwork  for  implementation 
of  the  USA-STAGS-CFA  system;  however,  some  additional  com¬ 
ments  are  required  as  well  as  a  reiteration  of  the  interconnection 
between  USA,  STAGS,  and  the  CFA.  With  regard  to  the  start¬ 
ing  procedure  outlined  in  §2.7,  it  should  be  noted  that  a  discon¬ 
tinuous  wavefront  cannot  be  propagated  “as  is”  through  the  fluid 
volume  mesh.  Rather,  the  wavefront  must  be  “ramped”  so  that 
its  value  at  the  front  is  one-half  of  the  jump,  in  line  with  the 
Fourier  convergence  theorem.  Only  with  this  modification  will  a 
discontinuous  wave  propagate  correctly. 
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Control  of  the  transient  analysis  is  governed  by  USA  with  STAGS 
and  the  CFA  functioning  as  subroutines  in  the  staggered  solution 
strategy: 

•  Convert  current  pressures  on  structure  boundary  to  forces  and 
obtain  structural  solution  with  STAGS.  Extrapolate  displacements 
to  f-f-  At.  If  the  structural  behavior  is  linear  the  trapezoidal  rule 
is  used  in  STAGS  and  the  extrapolation  is  that  of  (2.69).  If 
the  structural  behavior  is  nonlinear  the  Park  method  is  used 
in  STAGS  and  the  extrapolated  value  is  taken  to  be  the  last  value. 

•  Determine  displacements  on  DAA  boundary  due  to  incident  and 
hydrostatic  pressure  terms  at  t-j-At,  transform  current  pressures 
on  DAA  boundary  to  control  point  values  and  estimate  incremen¬ 
tal  displacements  due  to  total  pressure  at  t  +  At  from  (2.169). 
Sum  total  estimated  displacements  on  DAA  boundary. 

•  Using  as  input  the  estimated  displacements  on  the  structure  and 
DAA  boundaries,  solve  for  the  fluid- volume  pressures  at  t-j-  At 
using  the  CFA. 

•  Correct  the  DAA  boundary  displacements  using  the  new  total 
pressures. 

•  Save  system  responses  and  repeat  cycle. 
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§3.2  Usage 

Although  the  primary  emphasis  of  this  work  has  been  on  the 
development  of  the  CFA  and  its  interaction  with  USA-STAGS 
during  a  transient  response  analysis,  a  number  of  changes  have 
been  required  in  the  processing  that  must  preceed  the  time  in¬ 
tegration  phase  of  the  computations.  The  following  is  a  brief 
discusion  of  all  of  the  steps  involved. 

•  STAGS  Preprocessing.  Construct  the  structural  model  and 
create  grid  geometry  file.  Although  the  DAA  equations  used 
with  the  CFA  are  not  augmented  and  the  structural  mass  is  no 
longer  required  for  USA  preprocessing  the  structural  mass  file  is 
still  generated  because  it  also  contains  the  node-point/degree-of- 
freedom  information  necessary  for  USA  to  apply  pressure  forces 
to  the  structure  at  their  proper  locations. 

•  CFA  Preprocessing.  Construct  a  file  containing  the  following 
fluid  volume  information:  node  point  coordinates,  node  connec¬ 
tions  to  structure  and  DAA  boundaries,  node  constraint  tags,  and, 
an  element  node  list.  It  is  important  to  note  that  the  surface 
grid  on  the  structure  and  the  fluid  volume  grid  in  contact  with 
the  structure  should  be  identical  although  the  node  numbering 
schemes  need  not  be  the  same.  In  addition,  the  thickness  of  the 
volume  elements  in  contact  with  the  structure  must  be  carefully 
sized  to  meet  the  stability  criteria  developed  in  §2.5. 

•  USA  Preprocessing.  The  FLUMAS  processor  must  access  both 
the  structure  geometry  file  as  well  as  the  fluid  volume  file.  Al¬ 
though  fluid  control  points  must  be  defined  for  USA  on  the  DAA 
boundary  and  the  added-mass  matrix  created,  USA  must  keep 
track  of  the  fluid- volume/structure  connectivity.  Even  though 
augmentation  is  not  carried  out,  the  AUGMAT  processor  must  be 
executed  as  it  still  functions  to  produce  a  compact  data  base  file 
for  USA  to  access  during  the  transient  response  analysis. 
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•  USA-STAGS-CFA  Processing.  The  user  must  specify  a  level 
of  artificial  damping  for  the  CFA  and  a  “wavefront”  node,  i.e., 
the  fluid  volume  node  touched  by  the  wavefront  at  t  =  0,  which 
will  be  generally  be  located  as  close  as  possible  to  the  structure 
without  actually  being  on  it.  Although  the  incident  wave  could  be 
propagated  all  the  way  from  the  DAA  boundary  such  practice  is 
not  recommended  as  computer  time  is  simply  wasted  in  producing 
mesh  dispersion  that  erodes  the  wavefront  sharpness.  It  is  also 
important  that  the  values  chosen  for  the  artificial  damping  and 
time  stepsize  meet  the  stability  criteria  stated  in  §§2.5-2. 6. 

•  USA  Postprocessing.  Fluid  pressure  histories  can  be  obtained 
at  any  fluid  volume  node  desired  as  well  as  contour  plots  of  the 
pressure  field  at  specified  times.  A  consequence  of  this  capability 
is  an  increase  in  the  size  of  the  response  history  file  for  the  USA- 
STAGS-CFA  system  over  that  for  the  USA-STAGS  code. 
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Section  IV 


EXAMPLE  PROBLEMS 


§4.1  Overview 

The  USA-STAGS-CFA  system  has  been  tested  by  its  application 
to  two  problems  whose  solutions  have  been  obtained  by  other 
methods.  The  first  is  a  one-dimensional  problem  studied  by  Bleich 
and  Sandler  [5],  which  involves  a  flat  plate  initially  resting  on  the 
surface  of  a  half  space  of  fluid.  The  second  is  a  two-dimensional 
cylindrical  shell  surrounded  by  an  infinite  fluid  and  is  a  variant 
of  a  problem  discussed  by  Newton  [6].  For  both  problems  the 
excitation  consists  of  a  step-exponential  plane  wave  superimposed 
upon  an  ambient  hydrostatic  pressure  field. 

Bleirh-Sandler  P’ate  Problem 

This  is  effectively  a  one-dimensional  problem  whose  exact  solution 
can  be  obtained  by  the  method  of  characteristics.  The  USA- 
STAGS-CFA  model  consisted  of 

(1)  A  single  structural  square  plate  of  side  dimension  1.5  in.  and 
unit  thickness. 

(2)  100  cubical  fluid-volume  elements  of  side  dimension  1.5  in. 

(3)  A  DAA  boundary  with  a  single  control  point  at  the  center  of 
a  square  boundary  element  lying  on  a  face  of  the  last  volume 
element. 

Physical  properties  used  were  equivalent  to  those  of  [5];  however, 
they  were  converted  to  a  computationally  consistent  set  of  units. 
The  mass  density  of  the  plate  was  5.32986  x  10'4  lb  sec2  in-4 
while  that  of  the  fluid  was  9.3455  x  10~5  lb  sec2  in-4.  The  speed 
of  sound  in  the  fluid  was  57120  in  sec-1. 
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The  hydrostatic  pressure  in  the  fluid  at  the  plate  mass  was 
14.7  psi,  increasing  linearly  into  the  fluid  volume  consistent  with 
a  gravitational  acceleration  of  386.4  in  sec  2.  The  peak  pressure 
of  the  incident  wave  was  103  psi  and  its  decay  time  was  .9958  x 
10~4  sec. 

The  time  step  chosen  for  the  analysis  was  1.313  x  10-'5  sec 
(one  half  of  the  Courant  limit),  which  was  kept  constant  for  1200 
steps.  Four  sets  of  runs  were  made  with  and  without  cavitation 
allowed  and  using  artificial  damping  coefficient  values  of  (5  —■  0.0, 
0.25,  0.50  and  1.00.  According  to  the  stability  analysis  of  §2.5  the 
integration  process  should  be  unstable  for  =  0.0  but  stable  for 
the  other  three  values  and  this  was  in  fact  verified. 

Comparative  results  for  the  non-dimensional  upward  velocity 
of  the  plate  are  shown  for  the  stable  runs  in  Figures  6,  7,  and 
8.  Actual  velocities  in  in/sec  can  be  obtained  by  multiplying  by 
57.12,  while  the  time  scale  is  given  in  decay  time  units.  The 
solid  lines  are  the  USA-STAGS-CFA  results  whereas  the  discrete 
symbols  are  taken  from  the  solution  plots  in  [5j.  The  rapidly 
decaying  curves  that  are  essentially  zero  by  6  decay  times  are  for 
the  case  when  cavitation  is  not  allowed  while  those  that  continue 
out  to  12  and  beyond  illustrate  the  cutoff  effects  of  cavitation. 

As  can  be  seen  the  correlation  is  excellent.  The  results  show 
that  the  smoothing  effects  of  the  artificial  damping  have  only 
slight  influence  on  the  amplitude  and  timing  of  the  response.  The 
small  oscillations  perceptible  in  Figure  6  are  due  to  “ringing” 
and  dispersive  effects  of  the  fluid-volume  mesh,  and  gradually 
dissappear  as  is  increased. 
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Figure  8  Non-dimensional  upward  velocity  of  plate,  6 
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§4.3  Cylindrical  Shell  Problem 

In  [6]  Newton  presented  several  sets  of  results  for  an  infinite 
cylindrical  shell  problem  using  the  special-purpose  two-dimen¬ 
sional  code  DPLPOT,  but  a  direct  comparison  of  his  calcula¬ 
tions  with  those  of  USA-STAG S-CFA  could  not  be  readily  made. 
This  is  because  his  structural  model  included  a  damped  oscillator 
whose  properties  were  adjusted  to  simulate  a  neutrally  bouyant 
shell  with  internal  equipment.  The  oscillator,  however,  couples 
only  to  the  n  —  1  rigid-body  mode  of  the  shell.  As  USA-STAGS- 
CFA  is  not  based  upon  modal  superposition  and  could  not  easily 
emulate  the  oscillator  action,  Professor  Newton  kindly  consented 
to  rerun  one  of  his  cases  without  the  oscillator. 

The  structural  model  used  in  this  study  consisted  of  one  row 
of  twelve  STAGS  410-shell  elements  around  half  the  circumference 
in  order  to  take  advantage  of  symmetry.  The  shell  radius  was 
500  cm  and  the  axial  width  was  chosen  as  130.9  cm  so  that  the 
element  aspect  ratio  was  unity.  The  shell  wall  was  of  sandwich 
construction  with  2.5-cm-thick  face  sheets  separated  by  a  mass¬ 
less  core  29.1-cm-thick  that  was  allowed  to  carry  transverse  shear 
only.  The  physical  properties  of  the  shell  material  correspond 
to  structural  steel.  The  values  used  were  7.83  gm  cm-3  for  the 
density,  2.1  x  1012  gm  cm-1  sec'2  for  Young’s  modulus,  and  0.30 
for  Poisson’s  ratio.  The  behavior  of  the  structure  was  constrained 
to  be  linear  at  all  times. 

The  fluid-volume  model  consisted  of  192  brick  elements  ar¬ 
ranged  in  16  concentric  circular  cylindrical  layers  about  the  shell; 
each  layer  is  subdivided  into  twelve  equal  sectors  subtending  a 
15°  angle.  The  12-element  DAA  boundary  surrounds  the  fluid- 
volume  mesh  at  a  radius  of  2500  cm.  The  intermediate  radii  for 
each  sucessive  layer  increase  in  geometric  proportion  at  a  rate 
of  1.105823.  The  physical  properties  of  the  fluid  represent  sea 
water,  with  values  of  1.024  gm  cm-3  for  the  density  p  and  1.5  x 
105  cm  sec-1  for  the  speed  of  sound  c. 
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The  hydrostatic  pressure  was  107  gm  cm-1  sec-2  throughout 
the  fluid  volume.  The  peak  magnitude  of  the  step-exponential 
plane  wave  was  8x  107  gm  cm-1sec-2  while  the  decay  time 
was  5  milliseconds. 

The  time  step  used  in  the  transient  response  analysis  was 
.125  milliseconds  (the  Courant  limit  being  .352  milliseconds)  and 
the  artificial  damping  coefficient  p  was  taken  as  unity.  The 
response  calculation  was  carried  for  160  steps  out  to  a  time  of  20 
milliseconds  and  cavitation  was  allowed  to  occur  if  the  absolute 
pressure  dropped  below  a  value  of  zero. 

Comparative  results  are  shown  in  Figures  9  through  17.  The 
solid  lines  represent  the  USA-STAGS-CFA  calculations  while  the 
discrete  symbols  represent  Newton’s  computations  (which  were 
provided  at  half-millisecond  intervals).  Figures  9,  10,  and  11  show 
radial  displacement  histories  at  the  initial  point  of  contact  of  the 
wave  on  the  shell,  at  90°  from  the  contact  point,  and  at  180°  on 
the  back  of  the  shell,  respectively.  Figures  12,  13,  and  14  show 
radial  velocities  at  the  same  locations  while  Figures  15,  16,  and 
17  show  total  pressures  at  those  locations.  In  these  plots  the 
displacement  and  velocity  responses  have  been  rescaled  so  that 
the  length  measure  is  meters,  while  the  pressure  values  have  been 
multiplied  by  10~7  so  that  they  represent  megapascals.  Time 
is  given  in  milliseconds. 

As  can  be  seen  the  correspondence  is  very  good  for  the  selected 
displacement  and  pressure  histories,  except  for  the  displacement 
at  90°,  which  is  more  sensitive  than  those  at  0°  and  180°  to  dis¬ 
cretization  details  and  modal  convergence  problems.  Although 
the  USA-STAGS-CFA  velocity  responses  contain  oscillations  that 
do  not  appear  in  Newton’s  solution,  the  responses  agree  quite  well 
on  the  average,  especially  at  OP  and  180°.  (As  explained  in  the 
Remark  on  page  4-8,  these  oscillations  are  due  to  discretization 
effects,  and  do  not  have  physical  significance.) 

It  should  be  mentioned  that  the  early-time  pressure  peak  of 
Figure  15  is  a  real  effect  that  is  masked  in  the  DPLPOT  results 
by  the  coarse  output-sampling  interval. 
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It  is  remarkable  that  such  good  correlation  has  been  ob¬ 
tained  even  though  there  is  a  marked  disparity  in  the  discretiza¬ 
tion  details.  Newton  used  31  modes  for  the  structural  response, 
and  approximately  3000  two-dimensional  fluid-volume  elements 
of  similar  size  to  fill  a  rectangular  region  around  a  half  cylinder 
model.  Furthermore,  his  time  step  was  one-fourth  of  that  used 
in  the  USA-STAGS-CFA  computations. 

Although  the  structural  response  is  of  greatest  interest  in 
these  studies,  cavitation  does  occur  in  this  problem  and  an  idea 
of  its  extent  can  be  gathered  from  Figures  18  and  19,  which  are 
fluid-pressure  “snapshots”  at  8  milliseconds.  The  location  and 
shape  of  the  cavitating  region  is  roughly  the  same  in  both  sets  of 
computations;  the  region  closes  after  16  milliseconds. 

The  fact  that  structural  responses  agree  quite  well  despite  the 
use  of  a  much  coarser  fluid  model  in  the  USA-STAGS-CFA  com¬ 
putations  augurs  well  for  the  applicability  of  this  new  modeling 
capability  to  large  three-dimensional  underwater- shock  problems. 

Remark.  A  more  refined  USA-STAGS-CFA  analysis  of  this  problem  was 
carried  out  after  the  initial  draft  of  this  report  was  prepared.  The  spatial 
grid  was  halved  as  well  as  the  time  step,  resulting  in  a  320-step  calculation 
involving  768  fluid-volume  elements  and  24  structural  elements.  The  velocity 
oscillations  of  Figures  12  through  14  became  hardly  noticeable,  which  shows 
them  as  a  “coarse  grid"  effect.  The  displacement  and  velocity  histories  at 
$  =  QCP  displayed  better  correlation  with  Newton’s  results,  while  agreement 
for  the  other  sample  histories  remained  excellent. 
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Figure  18  Map  of  cavitation  zone  at  t=8  ms,  DPLPOT. 
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Figure  19  Map  of  cavitation  zone  at  t=8  ms,  USA-STAGS-CFA. 
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Code  69SG,  Y.  Shin 

Naval  Research  Laboratory 

ATTN: 

Code  8403,  R.  Bel  sham 

ATTN: 

Code  8440,  G.  O'Hara 

ATTN: 

Code  6380 

ATTN: 

Code  8100 

ATTN: 

Code  8301 

ATTN: 

Code  8406 

ATTN: 

Code  2627 

ATTN: 

Code  8445 

ATTN: 

Code  8404,  H.  Pusey 

Naval  Sea  Systems  Command 

ATTN: 

SEA- 08 

ATTN: 

SEA-55X1 

ATTN: 

SEA-033 

ATTN: 

SEA-06J,  R.  Lane 

ATTN: 

SEA-09G53 

ATTN: 

SEA-9931G 

ATTN: 

SEA- 323 

ATTN: 

SEA-0351 

Naval  Surface  Weapons  Center 

ATTN: 

Code  F34 

ATTN: 

Code  R13 

ATTN: 

Code  RIO 

ATTN: 

Code  U401,  M.  Kleinerman 
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Code  R14 

ATTN: 

Code  F31 

ATTN: 

Code  R15 

Naval  Surface  Weapons  Center 
ATTN:  W.  Wishard 
ATTN:  Tech  Library  t  Info  Svcs  B 

Naval  War  College 

ATTN:  Code  E-U 

Naval  Weapons  Center 

ATTN:  Code  343,  FKA6A2 
ATTN:  Code  266,  C.  Austin 
ATTN:  Code  3263,  J.  Bowen 

Naval  Weapons  Support  Center 

ATTN:  Code  70S53,  D.  Moore 
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Naval  Weapons  Evaluation  Facility 
ATTN:  G.  Binns 
ATTN:  Code  10 
ATTN:  Code  210 
ATTN:  R.  Hughes 

New  London  Laboratory 

ATTN:  Code  4494,  J.  Patel 
ATTN:  Code  4492,  J.  Kalinowski 

Newport  Laboratory 
ATTN:  Code  EM 


ATTN: 

Code  363,  P.  Paranzino 

Ofc  of  the  Deputy  Chief  of  Naval  Ons 

ATTN: 

OP  987 

ATTN: 

NOP  982,  Tac  Air  Srf  &  Ewdev 

ATTN: 

NOP  981 

ATTN: 

NOP  654,  Strat  Eval  &  Anal  B 

ATTN : 

OP  098T8 

ATTN: 

OP  982E,  M.  Lenzini 

ATTN: 

OP  957E 

ATTN: 

NOP  953,  Tac  Readiness  Div 

ATTN: 

OP  37 

ATTN: 

OP  225 

ATTN: 

OP  03EG 

ATTN: 

OP  21 

ATTN: 

NOP  952,  ASW  Div 

ATTN: 

OP  605D5 

ATTN: 

OP  981N1 

ATTN: 

OP  223 

Office  of  Naval  Research 

ATTN: 

Code  474 t  N.  Perrone 

Strategic  Systems  Project  Office 

ATTN: 

NSP-272 

ATTN: 

NSP-43 

ATTN: 

N5P-273 
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Air  Force  Institute  of  Technology 
ATTN:  Commander 
ATTN:  Library 

Air  Force  Systems  Conmand 
ATTN:  DLW 

*'r  Force  Weapons  Laboratory 
ATTN:  NTES-G,  S.  Melzer 
ATTN:  NTE,  M.  Plampndon 
ATTN:  NTES-C,  R.  Henny 
ATTN:  SUL 
ATTN:  NTED 

Assistant  Chief  of  Staff 

Intel  1 igence 
ATTN:  IN 

Ballistic  Missile  Office 
ATTN:  DEB 

Deputy  Chief  of  Staff 

Research,  Development,  &  Acq 
ATTN:  AFRDQ1 
ATTN:  R.  Steere 
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Deputy  Chief  of  Staff 

California 

Research  4  Technology,  Inc 

Logistics  4  Engineering 

ATTN 

S.  Schuster 

ATTN 

LEEE 

ATTN 

K.  Kreyenhagen 

ATTN 

M.  Rosenblatt 

Foreign  Technology  Division 

ATTN 

L  i  b  ra  ry 

ATTN 

NIIS  Library 

ATTN 

TQTD 

University  of  Denver 

ATTN 

SD8G 

ATTN 

Sec  Officer  for  J.  Wisotski 

ATTN 

SDBF,  S.  Spring 

Electric  Power  Research  Institute 

Rome  Air  Development  Center 

ATTN 

G.  Sliter 

ATTN 

RBES,  R.  Mair 

ATTN 

Commander 

Electro-Mech  Systems,  Inc 

ATTN 

TSLD 

ATTN 

R.  Shunk 

Strategic 

Air  Command 

General  Dynamics  Corp 

ATTN 

NRI-STINFO  Library 

ATTN 

J.  Mador 

ATTN 

J.  Miller 

OTHER  GOVERNMENT  AGENCIES 

ATTN 

M.  Pakstys 

Central  Intel  1 i gence  Agency 

Kaman  AviDyne 

ATTN 

OSWR/NED 

ATTN 

R.  Rue ten ik 

ATTN 

OSR/SE/F 

ATTN 

G.  Zartarian 

ATTN 

L  i  b  ra  ry 

NASA 

ATTN 

N.  Hobbs 

ATTN 

F.  Nichols 

ATTN 

R.  Jackson 

Kaman  Sciences  Corp 

ATTN 

Library 

DEPARTMENT 

OF  ENERGY  CONTRACTORS 

ATTN 

F.  Shelton 

University  of  California 

Kaman  Sciences  Corp 

Lawrence  Livermore  National  Lab 

ATTN 

D.  Sachs 

ATTN 

S.  Erickson 

Kaman  Tempc 

Los  Alamos 

National  Laboratory 

ATTN 

DAS  1  AC 

ATTN 

R.  Whitaker 

ATTN 

MS  S30,  G.  Spi 1 lman 

Karagozian 

and  Case 

ATTN 

Reports  Library 

ATTN 

J.  Karagozian 

ATTN 

MS  634,  T.  Dowler 

ATTN 

R.  Sanford 

Lockheed  Missiles  4  Space  Co,  Inc 

ATTN 

MS  670,  u.  Hopkins 

ATTN 

Technical  Information  Center 

ATTN 

T.  Geers 

Sandia  National  Lab 

ATTN 

B.  Almroth 

ATTN 
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Lockheed  Missiles  4  Spaces  Co,  Inc 
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DEPARTMENT 

OF  DEFENSE  CONTRACTORS 

M  4  T  Company 

ATTN 

D.  McNaight 

Applied  Research  Associates,  Inc 

ATTN 

D.  Piepenburg 

McDonnell 

ouglas  Corp 

ATTN: 

R.  Halprin 

Applied  Research  Associates,  Inc 

ATTN 

B.  Frank 

NKF  Engineering  Associates,  Inc 

ATTN 

R.  Belsheim 

BDM  Corp 

ATTN 

T.  Neighbors 

Pacific-Sierra  Research  Corp 

ATTN 

A.  Lavagnino 

ATTN: 

H.  Brode,  Chairman  SAGE 

ATTN 

Corporate  Library 

Pacifica  Technology 

Cal i fornia 

Institute  of  Technology 

ATTN 

A.  Kushner 

ATTN 

T.  Ahrens 

ATTN 

R.  Bjork 

ATTN 

G.  Kent 

Columbia  University 

ATTN 

H.  Bleich 

Physics  Applications,  Inc 

ATTN 

F.  Dimaggio 

ATTN: 

C.  Vincent 
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Physics  international  Co 
ATTN:  L.  Behrmann 
ATTN:  F.  Sauer 
ATTN:  J.  Thomsen 
ATTN:  E.  Moore 
ATTN:  Technical  Library 
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ATTN 

ATTN 

ATTN 

ATTN 
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T.  Cherry 
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D.  Grine 
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K.  Pyatt 
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Science  Applications,  Inc 

ATTN:  Technical  Library 

Southwest  Research  Institute 
ATTN:  A.  Wenzel 
ATTN:  W.  Baker 


SRI  International 
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W.  Wilkinson 
A.  Florence 


Westinghouse  Electric  Corp 
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Teledyne  Brown  Engineering 
ATTN:  J.  Ravenscraft 


Tetra  Tech,  Inc 

ATTN:  L.  Hwang 
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ATTN:  D.  Jortner 
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Weidlinger  Associates 
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